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ABSTRACT 



The complex LMS adaptive algorithm developed by Widrow 
[Ref. 1] is used in the frequency domain to estimate the 
azimuth and elevation angles of a plane wave incident upon 
a planar array. The complex LMS algorithm is applied to two 
cases. The first case is a passive detection problem. The 
second case is a pulse communication problem. In both 
cases, complex weights are determined using the complex LMS 
algorithm which cophase all of the output electrical signals 
from -the planar array. Three versions of the complex LMS 
algorithm are studied and their performances are compared. 
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INTRODUCTION 



I . 

Frequency domain beamforming is accomplished by applying 
appropriate phase shifts at the sensor outputs of an array 
to account for the relative propagation delays of a signal 
from a particular direction. The phase-shifted signals from 
all sensors are then added together coherently to realize 
the full array gain. Discrete Fourier Transform (DFT) beam- 
forming is the usual method of determining the direction of 
arrival of a plane wave signal. A discrete number of 
spatial frequency bins are formed and each bin corresponds 
to a discrete direction. If the number of spatial frequency 
bins is large, very fine spatial resolution can be obtained. 

The phase shifts needed to cancel the relative propaga- 
tion delays can be determined adaptively. The complex LMS 
adaptive algorithm is used in this thesis. The LMS adaptive 
filter adjusts its adaptive weights recursively to minimize 
the mean square difference between a reference signal and 
its estimate. When the beam is steered toward a signal 
propagating in a particular direction, the phase of the 
signals at all sensors must be the same. Therefore, the 
signal at any sensor can be used as a reference which the 
others will be matched. The estimated signal is obtained by 
weighting the input signal by the current adaptive weights. 
Note that no prior knowledge of the reference signal is 
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required. The response of the LMS adaptive filter converges 
to the discrete Wiener filter without a priori knowledge of 
the input [Refs. 1,2]. [Ref. 1] proposed the complex LMS 
algorithm to deal with complex inputs. [Ref. 5] addressed 
the implementation of the LMS adaptive filter in the fre- 
quency domain. [Ref. 4] used the LMS adaptive filter in the 
frequency domain to estimate the bearing of a plane wave 
due to an acoustic source radiating a sinusoidal signal. In 
this application, the LMS adaptive filter was implemented to 
estimate the phase difference between two sonar arrays 
separated by a distance many times the signal wavelength. 

The angle of arrival of a plane wave can be estimated if the 
frequency of the acoustic signal and the speed of wave propa- 
gation are known or can be extracted from the received signal 
itself. 

The objective of this thesis is to extend the results in 
[Ref. 3J and [Ref. 4] to a planar array of M xn elements 
(hydrophones) where M and N are greater than two. Such an 
array has an overall size many times the wavelength of the 
received signal. The inter-element separation, however, is 
usually maintained at a distance of less than or equal to 
one-half of the expected minimum wavelength. This requirement 
[Ref. 5J prevents the occurrence of grating lobes in the 
far-field beam pattern. A two-dimensional array allows 
spatial resolution in both azimuth and elevation. Even 
though the detection range in underwater acoustics is large 
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compared to the ocean depth, the effect of ray bending due 
to the inhomogeneous ocean medium can bend the incident 
acoustic wave such that it can appear to arrive at a steeper 
or shallower angle than the line of sight angle in the homo- 
geneous medium case. The elevation/depression angle is 
at present used to estimate Convergence Zone (CZ) and Bottom 
Bounce (BB) ranges. 

In this thesis, the problem of a plane wave incident upon 
a planar array of M xn elements is studied. The acoustic 
wave signal at each of the elements in the array are identical 
if the array is steered in the direction of the incident 
wavefront. However, if the main lobe of the array is not 
steered properly, the plane wave signal will have the same 
spectral content at each element but modified by a phase 
shift proportional to the location of the element with respect 
to some reference element. These undesirable phase shifts 
can be cancelled by applying appropriate phase weights at 
each element and thereby cophasing the total array output to 
realize its array gain. [Ref. 4] demonstrated that the LMS 
adaptive filter can achieve phase alignment between a refer- 
ence signal and input signal in the frequency domain by 
direct application of the complex LMS algorithm. This is 
equivalent to a tapped delay line structure in the time do- 
main. However, in the frequency domain, a time delay t 
corresponds to multiplication of a complex number that is 
equal to e^^^ where co is the signal frequency. The 
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implementation of the LMS adaptive filter in the frequency 
domain requires fewer computations per iteration than in the 
time domain. An added advantage of using phase weighting 
is that a continuous range of spatial directions can be 
described, whereas in a tapped delay line structure, only 
finite increments of delays can be applied. Figure 1 shows 
a functional diagram of an N-element adaptive array imple- 
mented in the frequency domain [Ref. 6] . 

Chapter II of this thesis describes the specific struc- 
ture of the adaptive filter and the equations implemented for 
simulations. The assumptions made in. the model are discussed 
and justified. Several modifications to the complex LMS 
adaptive algorithms are made to increase the array's spatial 
coverage and to ensure that the steady state phase weights 
do correspond to the direction of the incident wave. 

In Chapter III, a passive sonar system is modeled to 
test the ability of the modified complex LMS algorithm to 
estimate the direction of a source in the presence of noise. 

The simulation program is implemented in VS APL. 

Chapter IV demonstrates the application of the complex 
LMS algorithm to a pulse communication problem. The inte- 
gration time in this case is much shorter than that of the 
passive sonar case. Two types of pulse waveforms are included, 
continuous wave (CW) and linear frequency-modulated (LFM) . 

This simulation program is implemented in VS FORTRAN. 
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Sensor Pattern — forming network 




I J 



Figure 1. Frequency Domain Adaptive Array 
[Ref . 6 : p . 7 J 
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Chapter V concludes this thesis by identifying further 
research areas and other possible applications. 

Appendix A contains the derivation of the complex LMS 
algorithm. Appendix B has the description of the passive 
detection program implemented in APL. Appendix C has the 
description of the pulse communication program implemented 
in VS Fortran. 
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THEORY OF SYSTEM MODEL 



II . 

A. OVERVIEW OF ARRAYS [Refs. 5, 6,7,8] 

The characteristics of the array elements and their 
arrangement in forming the array determine the ultimate 
performance of an adaptive array system [Ref. 6]. Both the 
linear array and the planar array are examined here. 

1. Linear Arrays [Refs. 5,8] 

Consider a linear array that has M equally spaced, 
identical point source elements along the x-axis. For illus- 
tration, Figure 2 shows a 7-element linear array with uniform 
interelement spacing d^^ and a plane wave arriving at the 
array with an incident angle 9 as measured from the array 
normal. The phasor sum of all elements is: 



s(t) 



M-1 . , 

m=0 



( 2 . 1 ) 



where ip is the phase shift between the m^^ and the m+1^^ 
element for m = 0,1, 2, 3, ..., M-1 . 



ip = k 



d^ sin 9 



d sin 9 

2 TT / X \ 

C ^ 



( 2 . 2 ) 



where : 





Figure 2 . 



A Seven-Element Linear Array with Uniform 
Interelement Spacing and Point Source Elements 



k = wave number in radians per meter 

X = wavelength in meters 

f = frequency in hertz 

d„ = interelement spacing in meters 

A. 

9 = incident angle measured from array normal 

c = speed of wave propagation in meters per second. 



Rewriting equation (2.2) by substituting X = c/f yields 






27Tf 

c 



d^ sin 0 
c 



(2.3) 



The Fourier transform of equation (2.1) is: 



S(f) = F{s(t)} = / s(t)e ^^^^^dt (2.4) 

— oo 

00 M— 1 . • O jr 

= / I y (t)e^™'^e”^^'^^^dt (2.5) 

-00 m =0 

M— 1 . “ • -1 jr 

= / y (t)e"^^^^^dt (2.6) 

m=0 -oo 

s(f,0) = A(f,9)Y(f) (2.7) 



where Y(f) is the frequency spectrum of the incident wave and 
A(f,9) is called the space factor or array factor. The array 



19 



factor A(f,6) determines the directional plane of the array 
in a plane containing the array. The dependence of A(f,6) on 
frequency, speed of propagation, element separation, number 
of elements, and incident angle can be shown by rewriting 
A (f , 9 ) as : 

. ,27Tf , -ON 

M-l . , M-1 Dm(^d^sine) 

A(f,0) = I e^"'* = I e ^ (2.8) 

m=0 m=0 



Summing equation (2.8) yields: 

sin 

A(f,6) = e ^ (2.9) 

sin ^ 

The normalized directional pattern is given by: 

2 

G(f,6) = 10 log^Q{ (2.10) 

M 

For nonisotropic (non-point source) elements, it is necessary 
to introduce an additional factor E(f,0) in equation (2.9) 
to include the directional response pattern introduced by 
each sensor element [Ref. 6]. The overall directivity pat- 
tern then is given by the produce of the array factor- and the 
element factor [Ref. 5]. However, if the size of the individual 
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elements are small compared to a wavelength, they can be 
assumed to be omnidirectional point sources, i.e., E(f,0) = 1. 
The effects of increasing the number of elements while main- 
taining the same element spacing are shown in Figures 3 and 
4 [Ref. 6]. It can be seen that the main lob beamwidth de- 
creases as the number of elements increases, and the number 
of sidelobes and nulls increases. In Figures 5-8, the number 
of elements is held constant at M = 7 while the interelement 
spacing is varied to illustrate the effects of elements 
spacing on the directivity pattern. 

Beamsteering [Ref. 5] is accomplished by applying a 
linear phase shift across the line array as shown in Figure 
9 [Ref. 6] . The effect of the insertion of this sequence of 
phase shifts is that the main lobe is steered to an angle as 
measured off the boresight equal to 0 where 



0 



. -1 
sin 



1 
2 d 



6 } 



( 2 . 11 ) 



and 6 is the phase shift between adjacent elements. Figure 
10 shows the directivity pattern of a steered linear array. 

2. Planar Arrays [Refs. 5,6] 

Much of the analysis done in linear arrays can be 
extended to the case of a rectangular-shaped planar array. 

A circular planar array or a spherical volume array would re- 
quire the use of polar and spherical coordinates respectively. 
However, array theory is invariant under coordinate 
transformations . 
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Figure 3. 



-go*’ -45" 0" 45" 90" 

Azimuth angle 

Directivity Pattern of a Three-Element 
Linear Array [Ref. 6:p. 40 J 




-90" -45 O' 45' 90' 



Arimuth angle ** 

Figure 4 . Directivity Pattern of a Four-Element 
Linear Array [Ref. 6;p. 40] 
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Figure 5. Directivity Pattern of a Seven-Element 

Linear Array for d/A =0.1 [Ref. 6:p. 41] 



Azimuth angle 

- 90 * - 45 ° 0 “ 45 * 90 * 




Figure 6. Directivity Pattern of a Seven-Element 

Linear Array for d/A = 0.2 [Ref. 6:p. 41J 
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Figure 7. 



Directivity Pattern of a Seven-Element 
Linear Array for d/A = 0.5 [Ref. 6:p. 42J 
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Figure 8. 



Directivity Pattern of a Seven-Element 
Linear Array for d = A {Ref. 6:p. 42] 



Gam (d8) 




Array 

output 

y<t) 



Figure 9 . 



Beamsteer ing by Applying a Linear Phase 
Shift Across the Array [Ref. 6:p. 43] 
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Azimuth angle 




Figure 10. 



Steered Directivity Pattern iRef. 6;p. 44J 



A planar array has the advantage of resolving the 
azimuthal and the elevation angles of arrival of an incident 
plane wave [Refs. 5,9]. Consider a rectangular-shaped 
planar array as shown in Figure 11; sensor elements are 
arranged in a rectangular grid in the x-y plane. The center 
of the array is usually chosen as the coordinate origin. 

The entire array has M elements in the x-direction with uni- 
form spacing d„ and N elements in the y-direction with uni- 
form spacing d^. The elements are assumed to be point sources. 
The phasor sum of the entire array can be written as: 

■ jmip jmp 

“ II y(t) e ^ e ^ (2.12) 

m n 

where 

ip = 2 tt (—^) sind cos4> (2.13) 

X A 

dy 

ij; = 2 tt (^p) sinG sin({) (2.14) 

^y A 

The directivity pattern of the planar array is given by: 

jmip jmp 

A(f,e,4)) = [ [e y = A^(f ,0,(|))A (f ,0,(|)) (2.15) 

m n X y 

It follows from the model given by equation (2.12) that the 
planar array beam pattern is the product of the array factors 
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of two linear arrays. However, separability of the two- 
dimensional beam pattern is not necessary to ensure the 
proper operation of the LMS adaptive algorithm for a planar 
array. Beamsteering is accomplished by applying appropriate 
linear phase shifts for the row and column elements. For 
the rectangular array case under investigation here, it is 
more convenient to transform the elevation-azimuth (6,(j)) 
space to a rectilinear coordinate space (u,v) by the 
transformation : 



u = sin 6 cos (p (2.16) 

V = sin 0 sin (p (2.17) 

The parameters u and v are the direction cosines with 
respect to the x and y axes, respectively. Figure 12 shows 
the alternate diagrams for presenting two-dimensional array 
beam patterns [Ref. 6] . The ranges of the spherical angles 
are 0 £ 9 £ tt/ 2 and 0 ^ cp ^ 2 tt whereas the ranges of the 
rectilinear coordinate system are ~1 £ u £ 1 and “1 £ v £ 1. 

B. THE PLANE WAVE MODEL 

1 . Far-Field Condition 

The LMS adaptive filter designed in this thesis 
will provide spatial resolution for a planar array for an 
incident plane-wave field. The plane wave assumption is 
justifiable for a radiating source located in the far-field 
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Figure 11. Sensor Element Arrangement of a 

Rectangular Planar Array (Ref. 6:p. 45J 



^ V 








Figure 12. Transformation from Spherical Coordinate 

Space to Direction Cosine Space (Ref. 6:p. 46] 
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[Refs. 5,10], due to wavefront expansion. The far-field 
range for a planar array is given by [Ref. 5] : 



R > 




(2.18) 



where : 

L ^ L „ 1/2 

D = { + (^) } (2.19) 

represents the maximum radial extent of the transducer array 

and L and L are the dimensions of the planar array in 
X y 

the X and y directions, respectively. 

2 . Propagation of a Plane Wave from a Far-Field Source 
The plane wave solution of the Helmholtz wave equa- 
tion has the form: 



y(t,r) 



(27Tft+k*r) 



( 2 . 20 ) 



where y(t,r) is called the velocity potential, f is the 
frequency, k is the propagation vector, and r the position 
vector. In rectangular coordinates. 



k = kx + ky + kz 
— X y-^ z 



( 2 . 21 ) 



and 



/N /\ /S 

r = XX + yy + zz 



( 2 . 22 ) 
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For a planar array located at some reference location z = 0, 
the velocity potential is: 

j [27Tft+ (k x+k y) ] 

y(t,r) = y(t,x,y) = Ae ^ (2.24) 

For a planar array with discrete sensor elements located 

uniformly in the x-y plane with spacings d and d respec- 

tively, the continuous space variable x can be replaced by 

md„ and y replaced by nd . If the time signal is digitized 

for computer processing, then the time variable t can be 

replaced by ZT , where i is the discrete time index and T 

s s 

th.e sampling interval. To summarize: 

t ^ ilT (2.25a) 

s 

X ^ mdj^ (2.25b) 

y ^ ndy (2.25c) 

The corresponding discrete time, sampled space signal is 
given by: 



j2irfilT j (md k +nd k ) 

y(£T ,md^,nd ) = Ae ^e 

S X Y 



(2.26) 



where y(ilT ,md ,nd ) is usually shorted to y(il,m,n). 
S X Y 

propagation vector k at z = 0 has two components k 

— X 



The 
and k 

y 
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that can be related to the direction cosines u and v via 



[Ref. 5] ; 



and 



X 



2'im 

A 



(2 . 27a) 



2itv 

A 



(2 . 27b) 



Substituting equation (2.27) into equation (2.26) yields: 



y(£,m,n) 



j2-iTf£T j^(umd +vnd ) 
Ae A X Y 



(2.28) 



Let 



j 27Tf £T 

y(£) = Ae ® (2.29) 

represent the time dependence of the signal. Equation (2.28) 
then becomes : 

j^(umd +vnd ) 

y(£,m,n) = y(£)e (2.30) 

From equation (2.30), it is easy to see that the signal at 
each element location (m,n) has the same time dependence 
but has a different phase due to the element location and 
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the direction cosines associated with the incident angle of 
the plane wave upon the array. The exponential relationship 
of the phase in equation (2.30) also suggests that the 
phases in the x and y directions are separable. 

C. FREQUENCY DOMAIN LMS ADAPTIVE FILTER FOR SPATIAL RESOLUTION 
The objective of the LMS adaptive filter used in this 
thesis is to phase align- the signals from all sensor elements 
such that they add up coherently to realize the full array 
gain. Figure 13 shows a general cophasing scheme for linear 
arrays. Cophasing or phase alignment is done in the frequency 
domain by multiplying the frequency spectrum at each element 
by the proper phase weight in order to cancel out the phase 
due to element location. This is equivalent to a phasor 
rotation in the complex plane. The amount of rotation needed 
to align each sensor element is proportional to the frequency 
of interest, the direction cosines u, v and the location of 
that element. 

1 . Phase Weights for the Planar Array 

The total array output is maximized if all elements 
in the array are phase aligned. If we let cd(m,n) be the 
proper phase weight at location (m,n ) , then the phase 
weighted total array output is ; 

s(^) = II cd (m,n) y ( £ ,m,n) (2.31) 

m n 

Substituting equation (2.30) into equation (2.31) gives: 
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Figure 13. Adaptive Cophasing Scheme for a Linear Array 
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s(Z) = y(£) [ [ cd(m,n)e 




(2.32) 



m n 



If 



cd (m, n) 



-j — (umdj^+vndY) 



(2.33) 



e 



then the quantity inside the summations becomes unity and: 



where s(£) is the sum over all elements and equation (2.35) 
is the maximum signal level possible. This maximum level is 
achieved by tuning M xN adaptive weights cd(m,n) to conform 
to equation (2.33). The same phase weighting procedure is 
also true in the frequency domain, in fact, the implementation 
of phase weighting is inherently a frequency domain opera- 
tion. Since phase weight equation (2.33) is a function of 
wavelength X and X = c/f, the proper phase weight for co- 
phasing at each element is a function of frequency f. Thus 
for each valid frequency component in the signal, a different 



s(£) = y(£) I I (1) 



(2.34) 



m n 



or 



s(£) = MN y(£) 



(2.35) 
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set of phase weights {cd(m,n)} must be generated. Consider 
the DFT with respect to time of equation (2.32): 

j^(umd„+vnd ) 

S(q) = Y(q) I I cd(m,n)e (2.36) 

m n 

where q is the frequency index and Y (q) is the DFT of yil) . 

If equation (2.. 33) holds, then: 

S(q) = Y(q) I 1(1) = MN Y(q) (2.37) 

m n 

Only valid spectral lines will be processed. 

2 . The Frequency Domain LMS Adaptive Filter 

The general frequency domain LMS adaptive algorithm 
is derived in Appendix A. Suppose that the time signal z(Jl) 

/N 

is the reference or desired signal and z^(£) is the normalized 

sum of all signals in the planar array. In the frequency 

th 

domain, the reference signal in the q DFT bin is: 

L-1 

Z(q) = I zU)e ^ (2.38) 

Z=0 



and the estimated output in the frequency domain is : 



L-1 . -j^£q 

Z. (q) = I z.{Z)e ^ (2.39) 

5 ,= 0 ^ 
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^ 1 

Substituting equation (2.31) with z^(£) = yields: 



L-1 

Z^(q) = I ^ I I cd^(m,n)y(£,m,n)e (2.40) 

£=0 m n 

where : 

£ is the time index; £ = 0,1,..., L-1 
• q is the DFT bin index; q = 0,1,..., Q-1 

i is the complex phase weight iteration number. 

The DFT operation with respect to time in equation (2.40) 
can be performed first to yield: 

^i^^^ ^ ^ ^ cd . (m,n) Y (q,m,n) (2.41) 

m n 



where : 



Y (q,m,n) 



L-l 

I y(£,m,n)e 
£=0 



(2.42) 



Based on the complex LMS algorithm [Ref. 1], the adaptive 

t h 

filter output in the q bin is given by equation (2.41) . 

The error signal is generated by comparing the desired 
(reference) signal to the adaptive filter output (estimate) . 
The error is denoted by e^^, where: 
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e . 
1 



(2.43) 



= Z(q) - Z^(q) 

The estimate, equation (2.41), is formed by applying the 

t ll, 

phase weights {cd^(m,n)} of the i iteration to each element 
in the planar array. The complex weight cd^(m,n) is updated 
recursively as follows: 



cdi_|_i(m,n) = cd^(m,n) + 2y^e^y* (q,m,n) 



(2.44) 



where : 

m = 0,1,...,M-1 

n = 0,1,..., N-1 

(*) denotes complex conjugate 

y . = feedback coefficient, a parameter that con- 

trols the rate of convergence, algorithm 
noise, and the stability of the algorithm 
[Ref. 4] . 

t h 

From equation (2.44), it can be seen that the i+1 weight 
cd^_^j^(m,n) may have magnitudes larger than unity. This 
growth in magnitude is undesirable for the purpose of 
spatial resolution since spatial resolution depends on the 
relative phase between adjacent elements to resolve the 
direction of wave arrival. Thus, a normalization is neces- 
sary to bring equation (2.44) back to unity. This normali- 
zation is: 
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(2.45) 



cdi_^l(m,n) ^ (m,n)/ (m,n) 1 

These updated and normalized phase weights can now be applied 
to equations (2.41), (2.43), and (2.44) in sequence to 

compute the next set of phase weights. This iterative 
process stops when predetermined criteria are met. At that 
point, the set of phase weights {cd^(m,n)} can be used to 
find the direction cosines of the incident plane wave. 
However, the phase angles of the phase weights {cd^(m,n)} 
may have been wrapped around an integer multiple of 2ir. The 
procedure for phase unwrapping is explained in Sections 4 and 
5 of this chapter. 

The feedback coefficient in equation (2.44) con- 
trols the rate of convergence and the stability of the LMS 
adaptive filter. Robbins and Monro [Ref. 11] showed that 
the adaptive weights {cd^(m,n)} will converge to the optimum 

result if y . is allowed to decrease with the iteration index 
1 

i. The precise conditions are [Ref. 12] : 

y . >0 

1 

lim y . = 0 

1 

2_-^oo 




oo 
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< 



00 



I 

i=l 



y 



A coefficient y ■ that satisfies the above conditions will 

1 

work as long as the signal and noise inputs are truly 
stationary but will not be satisfactory for a filter operat- 
ing in a slowly varying environment. Widrow's LMS algorithm 
[Ref. 1] uses a constant value of y satisfying the inequality 



0 < y < 



A 



-1 

max 



where A is the largest eigenvalue of the correlation 
matrix of the input. Although this matrix is typically not 
known a priori, some bound can be set up by examining equation 
(2.44). If stationarity can be assumed, it is possible to 
update y^^ every iteration in order to obtain the optimum set 
of phase weights. In Appendix A a simple method of updating 
the feedback coefficient y^ is proposed to improve the per- 
formance of the LMS adaptive filter. 

3 . Applying the Frequency Domain LMS Adaptive Filter 
to a Planar Array 

Given a rectangular planar array with M xn elements, 
there are several ways to process the signal from the array 
to achieve spatial resolution in both azimuth and elevation. 
Three different ways are considered in this thesis. 
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a. Orthogonal Linear Arrays 

Two-dimensional spatial resolution is possible 
by just considering one linear array in the x-direction and 
one linear array in the y-direction. A total of M+N-1 ele- 
ments out of M X N are used. This scheme is useful when 
processing time is limited. Figure 14 illustrates the choice 
of the center linear arrays for this scheme. However, any 
two orthogonal linear arrays in the planar array can be 
used. The recursive equations needed to implement this 
algorithm are divided into two sets; one set for the linear 
array in the x-direction and the other for the y-direction. 

. In the x-direction the estimate is; 

. , M-1 

^x M ^ c (m) Y (q,m,n =n^) (2.46) 

i m=0 

where : 

n = constant y-direction index 
o ^ 

c^ (m) = unity magnitude phase weight. 

The error is the difference between the reference and the 
estimate . 



e 

X . 
1 



Z(q) - (q) 

i 



The recursive update for phase weights is: 



(2.47) 
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TV 




Figure 14 . A Particular Orthogonal Linear Arrays 
Configuration 
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(2.48) 



Ci^l (m) = c.(m) + Y*(q,m,n=n^) 

i i 



The phase of the new update is : 



Ci+i (iTi) I (2.49) 



In the y-direction, the procedure is similar: 



Estimate : 



(q) 



^ N-l 

I d^ (n) Y (g,m =m^,n) 
n=0 



(2.50) 



Error : 




Z(z) - Z (q) 
^i 



(2.51) 



Update: ‘^i+i “ d.(n) + 2y e Y*(q,m=m^,n) (2.52) 

1 1 ^i ^i ° 

Normalization: / 1 (2.53) 

The convergence constants and y^ are usually set to be 
equal since the statistics in the orthogonal directions of 
a planar array can be assumed to be the same for the obser- 
vation time of most systems. 

b. Two-dimensional Array 

This scheme uses all M x n elements and therefore 
it can realize the full array gain of equation (2.35). The 
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phase weights cd(m,n) are not assumed to be separable. The 
equations for the LMS adaptive filter in the frequency 
domain are: 

Estimate: Z^(q) = ^ I I cd^ (m, n) Y (q ,m, n) (2.54) 

m n 

/s 

Error: = Z(q) - Z^(q) (2.55) 

Update: cd^_^j^(m,n) = cd^(m,n) + 2y^e^Y* (q ,m, n) (2.56) 

Normalization: cd^^^(m,n) ^ (m,n)/ j (m,n) j (2.57) 

c. Separable Two-dimensional Array 

As mentioned in the discussion on planar arrays 
and plane waves, the form of a plane wave suggests that the 
phase of a signal at an element (m,n) is separable. This 
scheme then uses the separability property 

cd(m,n) = c(m)d(n) (2.58) 

to implement a two-dimensional LMS adaptive filter. All 
M xN elements are used but only M+N phase weights need to be 
updated recursively. The equations are: 

Estimate: Z^(q) = I d^(n) I c^ (m) Y (q,m,n) (2.59) 

n m 
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/\ 

Error; = Z (q) - Z^(q) (2.60) 

Updates: ^i+i ~ c^(m) + (n) Y (q,m,n) ) * (2.61) 

n 

di^l(n) = d.(n) + 2y^e^(J (m) Y (q,m,n) ) * (2.62) 

m 

Normalizations: (^) / I ^ I (2.63) 

*^i+l^^^ ^ (r^)/ (n) 1 (2.64) 



All three of the aforementioned schmes are implemented and 
their results compared. At the start of all three algorithms, 
the initial phase weights are set to the boresight of the 
planar array, i.e., magnitude equal to unity and phase angle 
to zero. The normalization of phase weights to unity forces 
the spatial transfer function of the LMS adaptive filter to 
have unit magnitude. The steady state phase response is 
designed to phase align all sensor elements in an element- 
by-element fashion. More discussion on this topic can be 
found in Appendix A. 

4 . Extracting Estimates of the Direction Cosines u and 
V from Phase Weights 

To extract u and v from the orthogonal linear arrays 
and the separable two-dimensional array cases discussed in 
Section C.3.a. and C.3.a.c., only M elements in the x-direc- 
tion and N elements in the y-direction need to be considered. 
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a. Direction Cosine Estimates for Linear Arrays 

Consider that in cases 3a. and 3c., the phase 

weights c (m) , d(n) have reached a steady state. The 

objective at this point is to relate the phase angles of 

these two sets of phase weights to their respective direction 

cosines. Let ^ (m) be the phase angle of c (m) and C (n) be 

X y 

the phase angle of d(n), i.e., 

/ N i- -1 r Im [c (m) ] -, 



and 



Cy(n) 



tan 



-1 



^Im [d (n) ] 
^Re[d(n)]^ 



(2 . 65b) 



It can be seen from equation (2.33) and using the concept 
of separability that: 



(m) = - •^(u md^) 



(2 . 66a) 



Cy(n) = - ^(v nd^) 

Solving equation (2.66) for u and v yields: 

^ -XC (m) 

u = — T ^ , m 7 ^ 0 

2Trmd^ ^ 

- -XC (n) 

= - 2irndy ' " ^ ° 



(2 . 66b) 



(2.67a) 



(2.67b) 
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Thus, in the x-direction, there are M-1 estimates of u; 

/s 

while in the y-direction, there are N-1 estimates of v. 

To find an estimate of the direction cosines from equation 
(2.67), one needs to take an arithmetic average of possible 
estimates ; 



u = 



1 Y 



M-l 



m=l 



27rmd 



(2.68a) 



X 



V = 



N-1 -Xe^(n) 

— z — ^ — 



N-1 



n=l 



27 Tnd, 



(2.68b) 



Equations (2.68a) and (2.68b) will be referred to as the 

/V /\ 

point by point method. Another way of finding u and v makes 
use of linear regression [Ref. 13]. Consider equation (2.66) 
where ^ (m) is a linear function of m with slope equal to 

27T^ 

— ^ linear function of n with slope equal 
2 TT ^ 

to — ^vd . Using a linear regression fit of M data points 

A Y 

vs. the element number m, the slope and intercept of the line 



^ (m) can be calculated. The same procedure can be used 

for 4^(n). Let the slope in the x-direction be s^^ and the 



slope in the y-direction be S 2 * 



Then 



2tt ^ , 

- — 



(2.69a) 



2 /\ 

- — ''■^y 



(2.69b) 
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/V /\ 

Thus, u and v can be solved by rearranging equation (2.69) 
to the form: 



u 



-As. 



27Td 



X 



(2.70a) 



V 



-As2 



(2.70b) 



The estimates u and v obtained from the linear regression 
method represent the best linear least-squares fit of the 
observed data. 

b. Direction Cosine Estimates for Planar Arrays 

The two-dimensional array discussed in Section 

3.b of this chapter does not require the phase weights to 

be separable. Consider that the set of phase weights 

{cd(m,n)} has reached a steady state. Let ^ (m,n) be the 

xy 

phase angle associated with the cd(m,n). Then: 



? (m,n) 

xy 



imicd (m,n) ] t 
^Re led (m,n ) ] ^ 



(2.71) 



From equation (2.33), it can be seen that 



*xy 



(m,n) 



- ^(u mdjj 



+ V 



nd^) 



(2.72) 



In general, 
and n. For 



5xy(ro,n) can be a more complicated function of m 
instance, in the near-field problem, one needs to 
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modify equation (2.72) to contain quadratic phase terms to 
account for wavefront curvature [Ref. 5]. For the plane 
wave model, equation (2.72) adequately describes the phase 
weights needed to steer the directivity pattern of the planar 

/V /N 

array to the direction corresponding to u and v. 

The point by point method is applicable here 
/\ /\ 

to find u and v given a steady state phase angle. However, 
a linear regression fit of ^^^(m,n) vs. the element index 
m and n appears to be more suited to this problem. Rewriting 
equation (2.72) in the form of the equation for a plane in 
three-dimensional space yields: 

C^^(m,n) = (^^^ud^)m + (-^vdY)n (2.73) 

2 TT ^ 

Equation (2.73) describes a plane with slope — ^ ud in the 

A X 

2 TT ^ 

m-direction (x-direction and slope — ^ vd in the n-direction 

A Y 

(y-direction) . Again, let the slopes be; 

^ u d^ (2.74a) 

®2 " ■ ^ ^ (2.74b) 

Thus, equation (2.74) is identical to equation (2.70) and 

/N /V 

the direction cosine estimates u and v are found by equation 
(2.71), i.e.. 
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(2.71a) 



u 



-Xs. 



2-nd 



X 



V 






(2.71b) 



This result is not surprising since the exponential represen- 
tation of the plane wave, equation (2.20), is inherently 
separable. 

5 . Unwrapping the Steering Phase Weights 
a. Linear Array Unwrapping 

The proper phase weights for beamsteering are 
given by equation (2.66) : 



C^(n» 



-27T 

A 



(u md^) 



(2 . 66a) 



5 (n) 
Y 



2 t \ , V 

^(v ndy) 



(2 . 66b) 



Consider a 7-element linear array lying in the x- 
with the center element as the reference element, 
index then rtins from m = -3, -2, -1, 0, 1, 2, 3. 
equal to 0.55 and d^ = A/2, then equation (2.67a) 



direction 

The element 
If u is 
reduces to: 



/V 

5 (m) = - TTum 



(2.75) 



= -0 . 55 7Tm 
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The following table shows the required phase weights ^ (m) 

/V 

needed to steer the beam to u = 0.55. 



TABLE 1 

PHASE WEIGHTS FOR BEAMSTEERING 



m -3 


-2 


-1 


0 1 


2 


3 


^ (m) 1.65 t7 


I.I77 


.5577 


0 - . 5577 


-I.I77 


-1.6577 


X 

?'(m) -.3577 


-.977 


. 55t7 


0 -.5577 


0.9t7 


.3577 


The phase factors 


{e 


} are a 


set of complex 


numbers 


in 



the complex plane. The angle of any complex number must lie 
within a 2 ir interval. The interval chosen here is [-17,77]. 
This means that any angle ^ (m) that is outside the range 
[-77,77] will be wrapped around to an angle (m) that is 
inside the range. This property can be shown as follows: 



^j(6+2T7k) ^ gjBgj 277k ^ 



since 



j 277 k 
e-* 



1 for k = 0,±1,±2,... 



therefore 



e 



jB 



has modulo 2 t7. 
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Referring back to Table 1, the observed angles C ' (m) are 
wrapped. If no processing is done to unwrap the observed 
angles, the equations derived in Section C.4 of this chapter 

/\ /V 

will not apply for all permissible values of u and v. For 

/V /N 

small values of u or v, the needed phase weights do not wrap 
around but the spatial range of interest is severely restricted 
[Ref. 4]. In tracking systems, the above restriction in look 
direction can be justified since a crude target direction is 
usually provided by a search array. The maximum spatial 
window of an M element linear array without the unwrapping 
of the steering phase weights can be calculated. For example, 
given- that the element spacing is d^^ = A/2, and the maximum 
permissible magnitude of ^ (m) is tt in the range [-tt,tt] , then 

equation (2.66a) becomes: 



> (m) 1 

’X 'max 






= I -TTum 



max 



(2.77) 



or 







1 



For the 7-element linear array described here, the maximum 
value that the index m can have is 3. Therefore, 



u 



max 



1 

3 
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In general. 




2 



(2.78) 



max 



M-1 



where M is the number of elements in the array. The corres- 
ponding angular coverage in terms of the polar coordinate 
angle 0 is (for ({) = 0) : 



The total angular cpverage is 20 . For the case of the 7- 
element linear array, this corresponds to a coverage of about 
40® out of a range of 180®. Figure 15 shows the expected 
phase angles required for beamsteering vs. element number. 
Figure 16 shows the wrapped angles vs. element number. It 
should be noted that either set of these angles (phase 
weights) will steer the beam to the proper direction. The 
difficulty with the wrapped (observed) angles is that the 

/N 

direction cosines u and v cannot be directly estimated using 
the methods developed in Section C.4 of this chapter. In 
order to unwrap the observed angles in Table 1, consider 
equation (2.76). It can be seen that the observed angles 
differ from the angles generated from equation (2.66) by 
an amount of ±27rk where k = 0,1,2,.... In order to unwrap 
the observed phase angles, it is necessary to find out which 
elements' phase angles have been wrapped around and by what 



0 



. -1 
sin 



(|u| 



(2.79) 



max 
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PHASE ANGLE IN RADIANS 




ELEMENT INDEX 



Figure 15. Phase Angles for Beamsteering 
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PHASE ANGLE IN RADIANS 




ELEMENT INDEX 



Figure 16 . Wrapped-Around Phase Angles 
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integer multiple of 2tt. 
that : 



Let w(m) be the unwrap factor such 



" ^x^^^ " ^x^^^ (2.80) 

/N 

Table 2 shows (for u = 0.55 and = X/2 ) how the implemen- 
tation of equation (2.80) will yield the phase weights 
required for spatial resolution. 



TABLE 2 

PHASE UNWRAPPING 



m 


-3 




-2 


1 


0 


1 


2 


3 


?^(m) 


1.657T 




1.177 


.5577 


0t7 


-.5577 


-1.177 


-1.6577 


5'(m) 


-.3577 




-.977 


.5577 


0 


-.5577 


0.9t7 


0.3577 


w(m) 


+2tt 




+2t7 


0 


0 


0 


-2T7 


-2t7 


?;(m) 


1.6577 




1.177 


.5577 


0 


-.5577 


-1.177 


-1. 65t7 


Observe from 


Table 2 and equation 


(2. 


66) that the phase 


angles for elements m = 


-1 and m 


= 1 


do not wrap around as 


long as d^ < 


A 

2 • 


Recall 


also that 


the 


center 


element 


is 


chosen as the 


: reference 


element . 


It 


is possible then after 



a set of steady phase weights has been computed that an 
estimate of the direction cosine can be computed using the 
phase angles at elements m = -1 and m = 1 in equation (2.67) 
below. 
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^ -ACy(m) 

u (m) = —7y — j — 

2iTmd 

X 



(2.67) 



For 



d 



X 



A/2, 



u (m) 



TT m 



(2 . 81) 



The estimate of direction cosine u using only the information 
(phase angles) at m = ±1 is denoted as u where 



= 4lu(l) + u(-l)] (2.82) 

g 2 



Equation (2.82) will yield a good estimate of u as long as 
_< The next step in this process of phase unwrapping 

is to use the result from equation (2.82) to generate a set 
of projected phase weights (m) } using the following 
relation : 



S^(m) 



-27Tmd^u 

X g 



A 



(2.83) 



Recall from equation (2.76) that a phase angle outside the 

rangle t-iT,Tr] is mapped to an angle within {-it, it] . By 

examining the magnitude of the projected angle, it is. possible 

to determine how many multiples of 2 it were lost due to the 

modulo 2 it property of complex numbers. The sign of u and 

y 
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the location of the element m determine the sign of the 
unwrap factor w(m). For linear arrays with an odd number of 
elements, the unwrap procedure can be described as follows: 

if (k-2)iT ^ I C (m) I < kiT (2.84) 

/\ 

and if | C (ni) | does not lie between the limits in equation 
(2.84), the value of k is decremented by two and the inequality, 

/V. 

equation (2.84), is tested again until ] ? (m) | falls within a 
2 tt interval, then 



1 w (m) I ■ = (k-1 ) 7T 



(2.85) 



For an M-element (odd) array, the initial value of k is 
M-1 

( 2 ) • The sign of w(m) is determined by: 



sign w(m) 



sign (m) 



sign (u ) 

g 



( 2 . 86 ) 



where : 

/M-1^ o n 1 /M-1^ 

m ( 2 / /•••/ 2 f i 

i.e., using the center element as the reference. A similar 
procedure can be utilized in an array with an even number of 
elements. The reference used should then be the point be- 
tween the two center elements since there is no need for the 
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reference point to coincide with the location of a sensor 
element. This choice of reference takes full advantage of 
the resulting symmetry. From Table 2 , it can be seen that 
only half of the elements need to be examined using equations 
( 2 . 84 ) - ( 2 . 86 ) . The unwrap factors {w(m)} for the other half 
are simply the negative of the first half. This unwrap 
procedure is good for the two array configurations discussed 
in Sections 3a. and 3c. of this chapter. The unwrap proce- 
dure for the two-dimensional array is similar but the compu- 
tation is a little more involved. The unwrapping procedure 
for the linear array in the y-direction is identical with 
the substitution d ^ d , m ^ n, ^ (m) ^ (n) , . . . , etc. 

X Y X y 



The unwrapping procedure is best illustrated by considering 



the following example. Consider the case (u,v) = (-0.7, +0.7) 
and d = d = A/2. Ecuation (2.87) then becomes: 

AY 



b. Two-dimensional Array Unwrapping 



The proper phase weights to steer a beam to 



(u,v) are given by: 





(2.87) 




-7T (urn + vn) 



( 2 . 88 ) 



-7T (-0 . 7m + 0 . 7n) 
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Assume also that we are given a planar array with M elements 
in the x-direction and N elements in the y-direction (the 
corresponding element indices are (x,y) (m,n) , with both 

M and N equal to 5. For the sake of symmetry, let the 
element indices run as follows; 



m 



M-1 M-1 M-1. 



n 



-( 



N-1 

2 



) ,-( 



N-1 

2 



)+l 



0 , 1 , . . . , (^) 



In this case, for M = 5, 



m = -2, -1,0, 1,2. 



Similarly, for the orthogonal direction. 



n = -2, -1,0, 1,2. 

The obvious choice of reference element is (m,n) = (0,0), 
i.e., the center element. The desired phase weights for 
this example are given in Table 3. 

The phase weights marked within the two triangles 
in Table 3 will be wrapped around, so the actual observed 
phases when the phase weights reach steady state are tabu- 
lated in Table 4 . 
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TABLE 3 



DESIRED PHASE WEIGHTS {5xy(rn,n)} GIVEN 
THAT u = -0.7,v = +0.7 



n 



2 


-2 . 8tt 


-2.17T 


• 

t — 1 


1 

o 

• 


0 


1 — 1 


-2.17T 


-1 . 4tt 


-0.7-it 


0 


0.7tt 


0 


-1 . 4it 


-0.7tt 


0 


0 . 7tt 


1.4tt 


1— 1 


-0.7tt 


0 


0.7tt 


1 . 4tt 


2 . Itt 


2 


0 


0.77T 


1 . 4tt 


2 . Itt 


2 . 8tt 




-2 


-1 


0 


1 


2 



TABLE 4 

OBSERVED PHASES OF THE STEADY STATE 
ADAPTIVE WEIGHTS (m,n) 

n 



2 


00 

• 

o 


-0 . Itt 


0 . 6tt 


-0 . 7tt 


0 


I — 1 


-0 . Itt 


0 . 6tt 


-0.7tt 


0 


0.7tt 


0 


0 . 6tt 


• 

o 

1 


0 


0.7tt 


-0 . 6tt 


r— 1 


-0.7tt 


0 


0.7tt 


-0 . 6tt 


0 . Itt 


2 


0 


0 .7tt 


-0 . 6tt 


O.lTT 


00 

• 

o 




-2 


-1 


0 


1 


2 • 



I 

I 
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Let w(m,n) be the two-dimensional unwrap factor, i.e., the 
unwrapped phase is : 



? (m,n) = (m,n) = (m,n) + w(m,n) 

^xy ' ^xy xy 



(2.89) 



By comparing (m,n) } 

xy 

4 , it can be seen that 
lated values in Table 5 



in Table 3 with (m,n)} in Table 

xy 

w(m,n) must be equal to the tabu- 



TABLE 5 

UNWRAP FACTORS w(m,n) 



n 

2 

1 

0 

-1 

-2 



-2tt 

-2tt 

-2tt 

0 

0 

-2 



-2tt -2tt 

-2tt 0 



0 0 
0 0 



0 2tt 



0 


0 


0 


0 


0 


2tt 


2tt 


2tt 


2tt 


2tt 



-1 



0 



m 



To generate w(m,n) , estimates of the direction cosines 



u and V must be computed using elements at m = ±1 and 



n = ±1. Let u and v be those estimates obtained using the 

g g 

point by point method below: 
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(2.90) 



u = j[u (1) + u (-1) ] 
g 2 

V = y[v(l) + v(-l)] (2.91) 

g z 



For d„ = d = \/2, equation (2.88) can be applied to com- 

A Y 

pute the projected phases (m,n)} by the relation: 

xy 



r (m,n) = - 7 t{u m + V n} for all m and n (2.92) 

^xy g g 



The set of projected phases (angles) are then examined to 
decide the proper unwrap factor for a particule element 
(m,n) . The logic is as follows: for each n, all elements m 

are examined. 



Check 



(k-2)TT < \E (m,n) < kir 

— ' ^xy ' ' 



(2.93) 



and if (m,n) [ does not lie between the limits, then k is 

xy 

decremented by two and equation (2.93) is tested again. If 



E (m,n) I does lie between the limits, then 
^xy ' ' i • 



w (m,n) I = (k-1) 7T 

xy 



(2.94) 



This procedure applies when M is odd and the initial value 
M-1 

of k is (-^) . The sign of w(m,n) is determined by: 

/N /N 

sign[w(m,n)] = -sign [sign (m) sing (u_) +sign (n) sign (v )] (2.95) 

y g 
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where : 



m 



,M-1 
^ 2 



. ,- 1 , 0 , 1 , 



.M-1 
^ 2 



) 



(N-1, T ^ T ,N-1, 

n — f m • • f "“-L ^ ) 



For the next value of n, equations (2.93), (2.94), and (2.95) 

are repeated for all values of m. This continues until the 
last value of n is reached. It is possible though to examine 
only half of the planar array since the unwrap factors for 
the other half are the negative of that of the first half. 

To ensure symmetry about the reference element, it is neces- 
sary to rotate the phase angle at the reference element to 
zero. This can be accomplished by multiplying the phase 
weights of all M N elements by the complex conjugate of 
the reference phase weight, i.e., 

cd . (m,n) cd . (m,n) cd . (m ,n ) (2.96 

where m and n are the indices locating the reference element. 
This operation will ensure that the unwrapping procedure 
will work properly. For linear array phase weights, this 
phase centering should also be completed before unwrapping. 
The equations are: 



c . (m) c . (m) c* (m ) 

1 1 1 o 



(2.97) 
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(2.98) 



d^(n) ^ d^(n) d^(n^) 

where is the index locating the reference element of the 
linear array along the x-direction and n^ is the index 
locating the reference element of the linear array along the 
y-direction . 

6 . Summary 

The original complex LMS adaptive filter [Ref. 1] 
requires three necessary modifications to make it useful for 
estimating the direction cosines of an incident plane wave. 
Without the following modifications, accurate spatial resolu- 
tion is not possible: 

normalization of the adaptive complex weights (phasors) 
to unity magnitude after each iteration. 

unwrapping of the observed steady state phase angles 
to extend the spatial coverage to the full range of 
u and V. 

allowing the feedback coefficient y to decrease with 
increasing iteration number to achieve convergence to 
an optimal set of phase weights and to realize a 
robust filter. 

D. NOISE MODEL 

In the SONAR environment, the ambient noise field is a 
composite of many different noise sources. Therefore, using 
the Central Limit Theorem [Ref. 13], the ambient noise can 
be modeled quite adequately as additive white Gaussian noise. 
Intentional jamming is not considered in this thesis. How- 
ever, the use of an adaptive filter to place a null at the 
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spatial location of a jamming signal has been studied exten- 
sively [Refs. 2,4,15] . The noise corrupted received signal 
is given by; 



r(£,m,n) = y(£,m,n) + n(£,m,n) 



(2.99) 



where y(£,m,n) is the sampled signal represented by equation 

(2.26) and n(£,m,n) is white Gaussian noise time samples 

2 

with zero mean and noise power If A is the signal 

amplitude (see equation (2.26)), then the signal-to-noise 
ratio is given by: 



SNR 




( 2 . 100 ) 



or 



(SNR)^g = 10 log (SNR) dB (.2.101) 

The performance of the frequency domain LMS adaptive filter 
was tested for various signal-to-noise ratios. The noise 
environment is assumed to be stationary during the look 
interval of the SONAR system. 
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III. LOW FREQUENCY PASSIVE SONAR TARGET LOCALIZATION 



The objectives in passive SONAR are to detect and possibly 
classify noise sources in the ocean. Most of the noise gener- 
ated by vessels is concentrated in the low frequency range 
(30-lK Hz). Acoustic energy in this frequency range is 
capable of long range propagation and, as a result, most 
long range detection systems in the underwater environment 
operate at these low frequencies. Propeller cavitation, 
machinery noise, and wake are the major sources of such noise. 
Under certain situations, very long range detection has been 
demonstrated in this frequency range. In a long range 
detection scneario, knowledge of the elevation/depression 
angle is necessary to resolve potential range ambiguities 
in convergence zone (CZ) problems. Therefore, the use of a 
planar array is well justified. 

A computer program was implemented using the VS APL 
language to test the performance of the LMS adaptive filter 
in a noisy passive environment. The array size used in the 
simulation is small compared to most modern systems but the 
other parameters are set to simulate a very realistic pas- 
sive SONAR. The APL language is used because of its large 
library of advanced signal processing functions and its 
interactive mode of operation which allows for rapid program 
development. The major disadvantage of the APL language is 
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its slower speed of execution compared to running a compiled 
Fortran program that performs the same functions. This is 
not a serious problem, however, for the research of this 
thesis . 

A. PROBLEM STATEMENT 

A low frequency acoustic signal from a far-field source 
is received by a planar array in a noisy underwater environ- 
ment. The noise is assumed to be Gaussian and uncorrelated 
with the signal. It is also assumed that the noise between 
elements is uncorrelated. The planar array has a square 
structure with M = 5 elements in the x-direction and N = 5 
elements in the y-direction. The element spacing is set to 
one-half the wavelength of the maximum frequency of the 
system's operating range. The geometry of the problem is 
shown in Figure 17 [Ref. 18] . The system parameters are: 

Signal: A cos(27Tft), where A is the amplitude and 

f is the frequency 

Integration time: 0.5 seconds (T^) 

- Frequency resolution: 2 Hz (1/T^ = f^) 

7 

Number of samples: 128 = 2 (L) 

- Number of sensor elements in the x-direction: 5 (M) 

Number of sensor elements in the y-direction: 5 (N) 

Speed of sound: 1500 meters/second (c^) 

Sampling rate: 256 samples/second (f ) 

s 

Frequency range: 0-128 Hz 

Number of frequency bins: 128 
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INCIDE 







Figure 17. System Geometry iRef. 18 :p. 19] 
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Elements spacing: 7.5 meters (d^ and d equal) 

A X 

2 

Noise model: additive white Gaussian ~ W(0,a^) 

Number of noise samples per sample function: 3200, 

i.e. , L xM X N = 128 x 5 x 5 

The signal used in the simulation is a 100 Hz sinusoid 
corrupted by white Gaussian noise. A total of 128 samples 
are taken for each processing period of 0.5 seconds. This 

5 

corresponds to 256 samples per second which satisfies the 
Nyquist sampling theorem [Refs. 12,16,17]. The maximum 
observed frequency in this case is 128 Hz. The 100 Hz 
signal will center in bin number 50 for this simulation. 
Frequency bin numbers 64-127 correspond to the negative 
frequencies [Ref. 16]. The element spacing of 7.5 meters 
is the maximum allowable separation for the 100 Hz signal 
to avoid grating lobes in the far-field directivity beam 
pattern. The speed of sound is the speed in the proximity 
of the planar array. 

B. SIMULATION 

Given the system parameters stated in Section A, and a 
128 point DFT , a 100 Hz sinusoidal signal will be centered in 
frequency bin number 50. The logical flow graph of the simu- 
lation program is shown in Figure 18. The principle of 
superposition allows different frequency bins to be processed 
independently. However, if the same frequency is emitted 
from two or more spatial locations, the LMS adaptive filter 
will lock on to the one closest to boresight. Complete 
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Figure 18. Logical Flow Graph of the Simulation 
Program for the Low Frequency Passive 
Detection Problem 
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documentation of the simulation program for this case can 
be found in Appendix B. The outputs of the simulation pro- 
gram are the estimated direction cosines u and v. The RMS 
error is defined as: 



£ 




1/2 

+ Av ) 



(3.1) 



where : 

/N 

Au = u - u 



(3.2) 



/N 

Av = V - V 



(3.3) 



where u and v are the actual direction cosines. This measure 
of error is consistent with the least-squares criterion used 
in estimating u and v. If the estimate of the spherical 

/N /S 

coordinates (6/(^)) are required, the following transformations 

/N /N /N /s 

will transform (u,v) to (6,4)) [Ref. 5]: 



6 (u, v) 



sin ^ [ (u) ^ + (v) 



1/2 



(3.4) 



and 



/N /N /N 



-1 



(J)(u,v) = tan (v/u) 



(3.5) 



It should be noted that the transformation from (u,v) to 

/N /N 

(6,(j)) is nonlinear. A particular value of the RMS error 
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in equation (3.1) can be due to infinitely many different 
values of 0 and 4> , that is, given an RMS error e, there is 
no unique set of spherical angle estimates. The spherical 
angle estimates (0,(1)) depend on the values of both u and v. 
Evaluating the total differential of equation (3.4) and 
equation (3.5) yields: 



d0 



/S /V. /\ /s. 

udu + vdv 



V 



^2 ^2 
(u+V^) 



(1 - 



^2 2 
(U+V ^) ) 



(3.6) 



and 



d({) 





(dv 



V du V 




(3.7) 



Replacing d0 by A0 , d({) by A(p , du by Au and dv by Av results 
in the following: 



A0 



uAu + vAv 



(u"^+v ) (1 - (u"^+v ) ) 



V 



A(j) 




u 



-2 
+ V 



(Av 



vAu . 

u 



(3.8) 



(3.9) 



C. RESULTS 

Three versions of the frequency domain LMS adaptive 
filter were discussed in Chapter II. The performance of 
each will be presented here. 
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1 . Orthogonal Linear Arrays 



This configuration uses only M+N-1 elements. The 
parameters used for the simulation are as follows : 

Signal - cos 2Tr(100)t (bin number 50) 

Number of iterations (I) - 75 

Initial feedback coefficient (y ) - 0.0005 

Initial feedback coefficient (Py) “ 0.0005 

Scale factor for y and y (a ) - 0.909 (see Appendix A) 

Input SNR at single array element, i.e., SNR at FFT 
input - 0 dB 

Incident angles (0,<})) - (55°, 35°) (see Figure 17) 

The corresponding direction cosines (u,v) - (0.67101,0.46985) 



Figures 19 and 20 show the convergence characteristics of 

the complex LMS adaptive filter vs. iteration number i. The 

solid horizontal lines are the true values of direction 

cosines u and v. The oscillations in the beginning of the 

adaptation are due to the large initial values of y and y . 

X y 

The feedback coefficients y and y are scaled down recur- 

X y 

sively by the scale factor via the rule described in Appendix 
A which is rewritten below: 






(3.10a) 



and 



75 



DIRECTION COSINE U 




ITERATION NUMBER 



Figure 19. Convergence Characteristics of the LMS 

Adaptive Filter using the Orthogonal Linear 
Arrays Configuration in Estimating the 
Direction Cosine u at an Input SNR of 0 dB 
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DIRECTION COSINE V 
0.42 0.44 0.46 0.48 




ITERATION NUMBER 



Figure 20. Convergence Characteristics of the LMS 

Adaptive Filter using the Orthogonal Linear 
Arrays Configuration in Estimating the 
Direction Cosine v at an Input SNR of 0 dB 
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Py(i-l)ag 



(3.10b) 



Uy(i) 



Convergence to the steady state occurred in fewer than 35 
iterations in both cases. Figure 21 shows the error in 
estimating u and v using equation (3.1) vs. iteration number. 
The use of a constant feedback coefficient requires a much 
longer iteration period and results in less accuracy (see 
Appendix A) . The chosen initial values of and y^ are 
outside the bound of the convergence coefficient described 
in Appendix A. However, the convergence coefficients are 
decreased rapidly using equation (3.10). This choice of 
y^ and y^ shows that the LMS adaptive filter with decreasing 
convergence coef f icient (s) is very robust. Table 6 summarizes 
the simulation results for a particular sample function of 
noise . 



2 . Two-dimensional Array 

This version uses all M xN elements and the complex 
weights are not assumed to be separable. The parameters 
used in the simulation are: 

Signal - cos 27r(100)t (bin number 50) 

Number of iterations (I) - 75 



Initial convergence coefficient - 0.001 

scale factor (a ) - 0.909 

s 

Input SNR at single array element, i.e., SNR at FFT 
input - 0 dB 
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RMS ERROR 




ITERATION NUMBER 



Figure 21. The RMS Error in Estimating u and v Using 
the Orthogonal Linear Arrays Configuration 
at an Input SNR of 0 dB 
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TABLE 6 



SUMMARY OF COMPLEX LMS ADAPTIVE FILTER 
PERFORMANCE FOR THE ORTHOGONAL LINEAR 
ARRAY CONFIGURATION AT SNR = 0 dB 



Simulation Values 





u 




0.671010 




V 




0.469846 




/\ 

u 




0.695675 




/\ 

V 




0.466038 


Au 


= u - u 




2.4665 xio”^ 


Av 


> 

1 

< > 
II 




-3.808 xio"^ 


e 


= (Au^+Av^) • 




0.024957 




6 




55® 




6 




54.783® 




4> 




35® 




/N 

4> 




33.818® 



Incident angle (6,4)) - (55®, 35®) 

Corresponding (u,v) - (0.67101,0.46985) 

Figures 22 and 23 show the convergence characteris- 
tics of the two-dimensional array configuration of the 
complex LMS adaptive filter vs. iteration number i. 

Figure 24 plots the RMS error in estimating u and 
V (see equation (3.1)). The results of the simulation for 
a particular sample function of noise are summarized in 
Table 7. 
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DIRECTION COSINE U 
0.60 0.65 0.70 0.75 




ITERATION NUMBER 



Figure 22. Convergence Characteristics of the LMS 
Adaptive Filter for the Two-dimensional 
Array in Estimating the Direction Cosine 
u at an Input SNR of 0 dB 
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DIRECTION COSINE V ' 

0.36 0.40 0.44 0.48 




ITERATION NUMBER 



Figure 23. Convergence Characteristics of the LMS 
Adaptive Filter for the Two-Dimensional 
Array in Estimating the Direction Cosine 
V at an Input SNR of 0 dB 
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RMS ERROR 

004 0.01 0.02 0.04 




Figure 24 . The RMS Error in Estimating u and v 

Using the Two-dimensional Array at an 
Input SNR of 0 dB 
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TABLE 7 



SUMMARY OF COMPLEX LMS ADAPTIVE FILTER PERFORMANCE 
FOR THE TWO-DIMENSIONAL ARRAY AT SNR = 0 dB 



Simulation Values 



u 


0.67101 


V 


0.46985 


/s 


u 


0.68442 




V 


0.46569 


Au 


1.341 xio"^ 


Av 


-4.16 xio"^ 


e 


1.404 xio"^ 


e 


55® 


/\ 


0 


55.876® 




35® 


/\ 




34.232® 



3 . Two-dimensional Array with Separable Weights , 

This configuration assumes that the complex weights 
are separable, i.e., cd(m,n) = c(m)d(n). The simulation 
parameters are the same as those used in the two-dimensional 
array case discussed in the previous section. Figures 25 
and 26 show the convergence characteristics of this configura- 
tion. Figure 27 shows the RMS error versus number of 
iterations. The summary for this run is in Table 8 . • An 
alternate way to illustrate the performance of the complex 
LMS adaptive filter is a plot of RMS error vs. input SNR at 



84 



DIRECTION COSINE U 




ITERATION NUMBER 



Figure 25. Convergence Characteristics of the LMS 
Adaptive Filter for the Two-dimensional 
Array with Separable Weights in Estimating 
the Direction Cosine u at an Input SNR of 
0 dB 
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DIRECTION COSINE V 
0.460 0.465 0.470 0.475 0.480 




ITERATION NUMBER 



Figure 26. Convergence Characteristics of the LMS 
Adaptive Filter for the Two-dimensional 
Array with Separable Weights in Estimating 
the Direction Cosine v at an Input SNR of 
0 dB 
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RMS ERROR 
0.001 0.01 




ITERATION NUMBER 



Figure 27. The RMS Error in Estimating the Direction 
Cosines u and v Using the Two-dimensional 
Array with Separable Weights at an Input 
SNR of 0 dB 
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TABLE 8 



SUMMARY OF COMPLEX LMS ADAPTIVE FILTER PERFORMANCE 
FOR THE TWO-DIMENSIONAL ARRAY WITH SEPARABLE 
WEIGHTS AT SNR = 0 dB 



u 

V 

/N 

u 

V 

Au 

Av 

£ 

0 

/s 

0 

c|> 

C|> 



simulation Values 
0.67101 
0.46985 
0.66304 
0.47353 
-7.97 xio"^ 
3.68 xio"^ 
8.78 xio"^ 

55® 

54.565® 

35® 

35.534® 



a single array element (SNR at FFT input) . Figure 28 shows 
these curves for all three configurations. 

D . SUMMARY 

The algorithm that assumes separable complex weights 
shows better performance over the prescribed SNR range. The 
improvement of performance of all three algorithms with 
increasing SNR is evident from Figure 28. Since each noise 
sample function has a total of 3200 independent samples, the 
RMS error versus SNR curve is rather smooth. If several sample 
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RMS ERROR 




SIGNAL-TO-NOISE LEVEL IN-DB 



Figure 28. The RMS Error in Estimating u and v 
Versus Input SNR. 

Curve 1 : Orthogonal Linear Arrays 

Conf iguration 

Curve 2 : Two-dimensional Array 
Curve 3 : Two-dimensional array with 
Separable Weights 
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functions are averaged, the slight 'humps' in Figure 28 
can be smoothed out further. The RMS error quantity is 
analogous to that of a 'miss distance' on a rectangular grid. 
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IV. PULSE COMMUNICATION [Ref. 18] 



The simulation results of two cases are presented in 
this chapter. They are: 

Homogeneous medium case, in which the output electri- 
cal signal data at each array element is produced by 
the ocean communication channel simulation computer 
’ program [Ref. 19]. Figure 29 shows the system geometry 
of this case [Ref. 18] . Note that the ray path from 
transmit to receive array is a straight path. 

Inhomogeneous medium case, in which the ray path is 
bent due to the variable sound-speed profile. Thus, 
the apparent direction of arrival viewed from the 
receive array is different from the previous case. 
Figure 30 shows the system geometry of the inhomo- 
-geneous medium case [Ref. 18] . 

From the analysis of Chapter II, it can be seen that the 
LMS adaptive filter should be able to phase align a planar 
array to point in the direction of arrival. The signals 
used for this simulation are a CW pulse and a LFM pulse. 
These are very common signals used in the SONAR environment. 
Information is carried by the modulation of these pulses. 

The Fortran program used to simulate this problem is docu- 
mented in Appendix C. Blount [Ref. 18] studied the effect 
of model-based cophasing on the probability of detection of 
a single pulse. The amount of cophasing is determined by 
the system geometry and deterministic ray bending. It was 
shown that by applying the phase weights generated by con- 
sidering those factors, the performance of a correlator 
receiver was improved markedly. Analysis done on those 
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Figure 29 . System Geometry for the Homogeneous 
Medium Case [Ref. 18:p. 16] 
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Figure 30. System Geometry for the Inhomogeneous 
Medium Case [Ref. 18 :p. 17] 
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steering phase weights showed that the main beam of the array 
directivity pattern was indeed steered to the direction of 
actual arrival instead of line of sight. 

A. TRANSMIT WAVEFORMS [Ref. 18] 

Two types of waveforms were used to test the LMS adap- 
tive algorithm. 

1 . Rectangular-envelope CW Pulse 

The signal presented to the processor is a quadra- 
ture demodulated complex envelope of the CW pulse [Refs. 5, 
18] . 

K ■ jk(27Tf )t 

z(t) = I ze ° (4.1) 

k=-K " 

where denotes complex envelope. 

The pulse repetition frequency is the same as the 
fundamental frequency f^ of the finite (K harmonics) fre- 
quency spectrum from which the pulse is ysnthesized. The 
pulse duty cycle is arbitrarily set to 0.5. The complex 
Fourier series coefficients z^ used to synthesize the com- 
plex envelope of the CW pulse are obtained from a closed-form 
expression for the complex-valued continuous spectrum. The 
Fourier coefficients are obtained by evaluating the closed- 
form expression for the continuous spectrum at discrete 
frequencies corresponding to integer multiples of the 
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fundamental frequency. The following specific transmit 
signal parameters were used in all CW pulse simulations: 
Amplitude (A): 40.0 

Duty Cycle (D) : 0.5 

Fundamental Frequency (f^) ^ 200 Hz 

Harmonic Values: = 20 exp [jO®] 

z_j^ “ ^1 “ 12.3324 exp [j0°] 
z _2 = Z 2 = 0.000 exp [jO®] 
z _2 = z^ = 4.244134 exp {jl80°] 
z_^ = z^ = 0.000 exp [j0°] 
z_^ = = 2.546479 exp [j0°] 

2. Rectangular-envelope LFM Pulse [Ref. 18] 

The complex Fourier coefficients used to synthesize 
the LFM pulse are found using a procedure similar to that 
used for the CW pulse except the closed form expression for 
the complex-valued continuous spectrum of the LFM pulse was 
found by using the method of stationary phase. Officer 
[Ref. 20] describes the method of stationary phase as does 
Papoulis [Ref. 21] who also provides a complete description 
of the LFM waveform. The following transmit signal param- 
eters were used in all LFM pulse simulations: 

Amplitude (A): 40.0 

- Duty Cycle (D) : 0.8 

Phase Deviation Constant (B) : 2356.2 radians/volt 

Fundamental Frequency (f ) : 10 Hz 
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Harmonic Values: = 14.60593 exp [j45®] 

z_^ = = 14.60593 exp [j21®J 

z _2 - z^ - 14.60593 exp [309®] 

z _2 = z^ = 14.60593 exp [jl89®] 

B. PROBLEM STATEMENT 

A pulsed signal (CW or LFM) is sent to an intended 
receiver in the far-field. The signal at the receive array 
has a planar wavefront. The direction of arrival of the 
incident plane wave is determined by applying the frequency 
domain LMS adaptive filter to phase align (cophase) signals 
at all sensor elements in the planar array. The discrete 
time signal at element (m,n) has the form; 

2 TT 

j- 5 -(umd„+vnd ) 

y(£T ,md ,nd ) = z(£T )e (4.2) 

Sax S 

where 



K jk2TTf t 

z(£T ) = \ z,e ° (4.3) 

® k=-K 

£ : time index 

m: element index in the x-direction 

n: element index in the y-direction 

T: sampling period 

u: direction cosine with respect to the x-axis 
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V 



direction cosine with respect to the y-axis 
d: interelement spacing in the x-direction 

d: interelement spacing in the y-direction 

z: complex Fourier coefficient 

f; fundamental frequency 
K: total number of harmonics 

[Ref. 5] shows that the number of time samples needed to 
completely describe equation (4.3) is L, where: 



L > 2K + 1 



(4.4) 



The system parameters are: 



CW 



Integration time (T^) : 

Frequency resolution (f^) : 

Number of samples (L) : 

Number of sensors (M) : 

Number of sensors (N) : 

Sampling rate (f ) : 

s 



5 mS 
200 Hz 
11 
5 
5 

2200 samples 
per second 



LFM 
100 mS 
10 Hz 
7 
5 
5 

70 samples 
per second 



Number of frequency 
bins (Q) : 

Element spacing d„, d^^: 

A Y 

Carrier frequency (f^) : 
Noise model: 



11 

0.1229 meters 
5 KHz 



7 

0.1229 meters 
5 KHz 



Additive white Gaussian noise for 
both cases 



Number of complex noise 

samples/pulse: 275 



175 
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The carrier frequency f is assumed to be known a priori. It 
is found for the CW pulse that the best direction cosine 
estimates are obtained by processing the spectral line corres- 
ponding to the carrier frequency. This is reasonable since 
the signal-to-noise ratio at that bin is the highest. For 
the LFM pulse, all the harmonic lines still have roughly the 
same magnitude after propagating through the medium. The 
accuracy in estimating the direction cosines is about the 
same for any harmonic line. Therefore the spectral line 
corresponding to the carrier is processed. 

C. RESULTS 

The homogeneous medium case is considered first for both 
CW and LFM pulses, followed by the inhomogeneous case. 

1 . Homogeneous Case 

The parameters for the system geometry (Figure 29) 

are : 

Speed of sound (c^) : 1475 meters/second 

Depth of transmit array (Y^) : 1000 meters 

Depth of receive array (Y^) : 2500 meters 

Cross range (X^ meters 

Line of sight range |r - r^ j : 3000 meters 

True spherical angle 9: 31.81° 

True spherical angle (p : -108.4° 

Direction cosine u: -0.1666 

Direction cosine v: -0.5000 
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The three array configurations of the frequency domain LMS 
adaptive filter were then applied to the CW and LFM pulse 
cases to estimate the direction cosines. The averaged re- 
sults for one pulse corrupted by 100 sample functions of 
noise at an input SNR of 0 dB at a single array element for 
the CW pulse, homogeneous case are presented in Table 9. 

The initial convergence coefficient is 0.5 and the scale 

factor a is 0.909 for both CW and LFM waveforms, 
s 

TABLE 9 

PERFORMANCE OF COMPLEX LMS ADAPTIVE FILTER FOR 
SPATIAL RESOLUTION, 100 ITERATIONS, INPUT 
- SNR = 0 dB FOR CW PULSE, HOMOGENEOUS CASE 

Algorithm 

orthogonal 2-dimensional 

linear arrays with separable 2-dimensional 







phase weights 


array 


u 


-0.1666 


-0.1666 


-0.1666 


V 


-0.5000 


-0.5000 


-0.5000 


u 


-0.1714 


-0.1672 


-0.1670 


V 


-0.4757 


-0.5003 


-0.5032 


Au 


-0.00048 


-0.0006 


-0.0004 


Av 


0.0243 


-0.0003 


-0.0032 


e 


0.0247 


-0.000624 


0.0032 


6 


O 
1 — 1 
00 

• 

1 — 1 
ro 


31.81° 


31.81° 


4) 


-108.4° 


-108.4° 


-108.4° 


6 


30.37° 


31.84° 


32.02° 


/N 

4> 


-109.8° 


-108.5® 


-108.4° 
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The rms error versus signal-to-noise ratio curves for all 
three array configurations are presented in Figure 31. 

The simulation results for the LFM pulse in the 
homogeneous medium at an input SNR of 0 dB at a single array 
element are summarized in Table 10. A total of 100 differ- 
ent noise sample functions were used to corrupt the received 
signal. The tabulated . values in Table 10 are the ensemble 
average values. Figure 32 illustrates the rms error versus 
SNR plot for the LFM pulse in the hom.ogeneous medium. 

TABLE 10 

PERFORMANCE OF COMPLEX LMS ADAPTIVE FILTER FOR 
SPATIAL RESOLUTION, 100 ITERATIONS, INPUT 
SNR = 0 dB FOR LFM PULSE, HOMOGENEOUS CASE 

Algorithm 

orthogonal 2-dimensional 2-dimensional 





linear arrays 


with separable 
phase weights 


array 


u 


-0.1666 


-0.1666 


-0.1666 


V 


-0.5000 


-0.5000 


-0.5000 


u 


-0.1096 


-0.1611 


-0.1462 


V 


-0.2695 


• -0.4096 


-0.4079 


Au 


0.0570 


0.0055 


0.0204 


Av 


0.2305 


0.0904 


0.0921 


e 


0.2375 


0.0905 


0.0943 


e 


31.81° 


31 . 81° 


O 
1 — 1 
00 

• 

1 — 1 
00 




-108.4° 


-108.4° 


-108.4° 


e 


16.91° 


0 

1 — 1 
1 — 1 

• 

(N 


25.68° 


/\ 


-112.1° 


-111.5° 


-109.7° 
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RMS ERROR 




SIGNAL TO NOISE LEVEL IN DB 



Figure 31. 



The RMS Error in Estimating u and v Versus 
Input SNR for the CW Pulse/Homogeneous 
Medium Case. 

Curve 1: Orthogonal Linear Arrays 

Configuration 

Curve 2: Two-dimensional Array 

Curve 3: Two-dimensional Array with 

Separable Weights 
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RMS ERROR , 




SIGNAL TO NOISE LEVEL IN DB 



Figure 32. 



The RMS Error in Estimating u and v Versus 
Input SNR for the LFM PUlse/Homogeneous 
Medium Case. 

Curve 1 ; Orthogonal Linear Arrays 
Configuration 

Curve 2: Two-dimensional Array 

Curve 3 : Two-dimensional Array with 

Separable Weights 
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The performance of the LFM pulse is worse than that 
of the CW pulse since the number of samples per pulse is 
seven where the CW pulse has eleven time samples per pulse. 
Consider also the Fourier coefficients given in Section A of 
this chapter. The magnitude of the spectral line correspond- 
ing to the carrier for the CW pulse is 1.37 times larger than 
that of the corresponding spectral line in the LFM pulse. 

However, for both waveforms, it can be seen from Figures 

31 and 32 that the two-dimensional array with separable 

weights has the best performance, followed by the two-dimensional 

array and orthogonal linear arrays, respectively. For the 

CW pulse, accurate spatial localization can be obtained for 

SNR greater than -3 dB . For the LFM pulse, the SNR required 

is about 3 dB. 

2. Inhomogeneous Case [Ref. 18J 

The parameters for the system geometry (Figure 32) 

are : 

Speed of sound (c^) : 1475 meters per second 

Gradient: 0.017 per second 

Depth of transmit array (Y^) : 1000 meters 

- Depth of receive array (Y^) : 2500 meters 

- Cross range (X^ meters 

Line of sight range |r - r^ | : 3000 meters 

- True spherical angle 0: 30.10° 

True spherical angle (|) : -109.4° 

- Direction cosine u: -0.1666 

- Direction cosine v: -0.4731 
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Figure 30 shows that the path traveled by the acoustic rays 
is bent. The LMS adaptive filter is able to resolve the 
actual direction of arrival but not the true line of sight 
direction to the transmit array. The initial convergence 
coefficient y is set equal to 0.5 and the scale factor a 

s 

is set equal to 0.909. Both the CW and LFM pulses are con- 
sidered for each of the three array configurations. Table 11 
summarizes the performance of the LMS adaptive filter applied 
to the CW pulse case at an input SNR of 0 dB at a single 
array element. Figure 33 shows the decline of rms error as 
the signal- to-noise ratio is increased. All three array 
configurations are included in the plot for comparison. 

The performance for the LFM pulse case is tabulated 
in Table 12 for an input SNR of 0 dB at a single array ele- 
ment. Figure 34 illustrates the performance of all three 
array configurations versus signal-to-noise ratio for the 
LFM pulse case. 

D. SUMMARY 

For both the homogeneous and inhomogeneous medium cases, 
the complex LMS adaptive filter performed as expected. The 
array configuration that assumed separability of the complex 
weights consistently demonstrated better performance than 
that of the orthogonal linear arrays and the two-dimensional 
array. The superior performance can be attributed to the 
fact that equation (2.24) which describes the reception of 
a plane wave by a planar array is separable. Therefore, by 
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TABLE 11 



PERFORMANCE OF COMPLEX LMS ADAPTIVE FILTER FOR 
SPATIAL RESOLUTION, 100 ITERATIONS, INPUT 
SNR = 0 dB FOR CW PULSE, INHOMOGENEOUS CASE 



Algorithm 





orthogonal 
linear arrays 


2-dimensional 
with separable 
phase weights 


2-dimensional 

array 


u 


-0.1666 


-0.1666 


-0.1666 


V 


-0.4731 


-0.4731 


-0.4731 


u 


-0.1663 


-0.1686 


-0.1651 


/N 

V 


-0.4527 


-0.4782 


-0.4746 


Au 


0.0003 


-0.0019 


0.0015 


Av 


0.0204 


-0.0051 


-0.0015 


z 


0.0204 


0.0054 


0.0022 


e 


30.10° 


30.10° 


30.10° 




-109.4° 


-109.4° 


-109.4° 


/s 

e 


28.83° 


30.47° 


30.17° 


4) 


-110.2° 


-109.4° 


-109.2° 
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RMS ERROR 




SIGNAL TO NOISE LEVEL IN DB 



Figure 33. The RMS Error in Estimating u and v Versus 
Input SNR for the CW Pulse/Inhomogeneous 
Medium Case. 

Curve 1: Orthogonal Linear Arrays 

Configuration 

Curve 2: Two-dimensional Array 
Curve 3: Two-dimensional Array with 
Separable Weights 
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TABLE 12 



PERFORMANCE OF COMPLEX LMS ADAPTIVE FILTER FOR 
SPATIAL RESOLUTION, 100 ITERATIONS, INPUT 
SNR = 0 dB, LFM PULSE, INHOMOGENEOUS CASE 



Algorithm 





orthogonal 
linear arrays 


2-dimensional 
with separable 
phase weights 


2-dimensional 

array 


u 


-0.1666 


-0.1666 


-0.1666 


V 


-0.4731 


-0.4731 


-0.4731 


/\ 

u 


-0.1035 


-0.1614 


-0.1417 


V 


-0.2447 


-0.4058 


-0.3332. 


Au 


0.0631 


0.0052 


0.0249 


Av 


0.2284 


0.0673 


0.1399 


e 


0.2369 


0.0675 


0.1421 


0 


30.10*^ 


O 

o 
1 — 1 

• 

o 

ro 


30.10° 


< 1 ^ 


-109.4° 


-109.4° 


-109.4° 


/s 

e 


15.41° 


25.89° 


21.23° 


/N 

< 1 ^ 


-112.93° 


-111.7° 


-113.0° 
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RMS ERROR 




SIGNAL TO NOISE LEVEL IN DB 



Figure 34 . 



The RMS Error in Estimating u and v Versus 
Input SNR for the LFM Pulse/Inhomogeneous 
Medium Case. 

Curve 1 : Orthogonal Linear Arrays 

Configuration 

Curve 2: Two-dimensional Array 

Curve 3 : Two-dimensional Array with 

Separable Weights 
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assuming separable complex weights, the system is set up to 
match the physics of the problem. The performance of this 
pulse communication system can be enhanced by increasing the 
pulse width, taking more time samples per unit time, using 
high resolution spectrum analysis, and enlarging the size 
of the array by adding more sensor elements. 
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V. CONCLUSIONS AND RECOMMENDATIONS 



The frequency domain LMS adaptive algorithm has been 
shown to perform the function of spatial resolution. The 
number of iterations required (approximately 35) to reach a 
steady state is found to be much fewer than that of a com- 
parable time domain adaptive filter [Refs. 2,4]. The three 
modifications made to the original complex LMS adaptive 
filter [Refs. 1,4] are: 

- normalization of the complex weights to unity 
magnitude after each update 
*• • 

reduction of the magnitude of the convergence 
coefficient for each iteration 

unwrapping of the phase weights 

It has been shown that the above three modifications enable 

the frequency domain LMS adaptive filter to be applied to 

a multiple element array, to have a fast convergence rate 

and robustness, and to cover the entire angular range of 

9 and (p values . 

In Chapter III, a passive low-frequency signal was 
generated to test the performance of the frequency domain 
LMS adaptive filter in the presence of white Gaussian noise. 
The low frequency problem is significant in the area of long 
range detection and possibly long range control and communi- 
cation for ballistic missile submarines. Some form of 
modulation, whether amplitude or angle, will have to be used 
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to carry the information. The slow data rate is not a major 
drawback since only certain predefined emergency operation 
codes will be transmitted. The simulation results indicate 
that for an input SNR at a single array element of greater 
than -3 dB, accurate bearing and depressive angle estimates 
can be obtained. 

In Chapter IV, a pulse communication problem was con- 
sidered. Two types of pulses, rectangular-envelope CW and 
rectangular-envelope LFM pulses, were used. The possible 
applications at this frequency range (5 KHz) are high reso- 
lution SONAR imaging, active SONAR detection, and communica- 
tion -between submerged vessels. The performance of the 
frequency domain LMS adaptive filter was again tested in the 
presence of white Gaussian noise. The simulation results 
indicate that accurate bearing and depression angle esti- 
mates can be obtained for input SNRs of greater than 0 dB . 
The result here is slightly worse than that of the passive 
low-frequency case due to the shorter data length used for 
the pulse communication waveforms. For the simulations 
in Chapter IV, the number of time samples collected is the 
minimum required to satisfy the Nyquist sampling theorem. 

The following three configurations of the frequency 
domain LMS adaptive filter were considered: 
orthogonal linear arrays 

- two-dimensional array with separable phase weights 
two-dimensional array 
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In all simulations, the configuration that used the separa- 
bility assumption produced the best results. This is quite 
reasonable since this configuration uses all the output data 
from the available elements and more importantly, its assump- 
tion of separability matches the physics of the plane wave 
signal since the phases of a plane wave in the orthogonal x 
and y directions are separable. 

In the course of this investigation, some interesting 

topics for further research revealed themselves: 

exact target localization (range, depth, bearing, 
and depression angle) if accurate environmental 
data can be obtained 

- -more efficient update algorithm for the convergence 
coefficient 

implementation of more efficient software 

application of the steady state phase weights generated 
by the frequency domain LMS adaptive filter to Blount's 
[Ref. 18] correlator receiver 

addition of a noise reduction system before spectral 
estimation 

- modifications to make the system jam-resistant 

applications of other high resolution spectral analysis 
techniques to produce frequency spectra (such as 
maximum entropy [Ref. 22], autoregressive, maximum 
likelihood, etc.). 
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APPENDIX A 



THE FREQUENCY DOMAIN LMS ADAPTIVE ALGORITHM 
A. DERIVATION 

Consider the phase aligning system shown in Figure 13. 

The performance objective is for the adaptive filter to con- 
verge to a set of phase weights such that the array output 
signal will match a reference signal. In the frequency 
domain, the output of a linear array of M elements at iter- 
ation i can be written as: 

2i(q) = C^(m)Y(q,m) = ^ c^Y(q) (A.l) 

m 

If Z (q) is the reference signal, then the error between the 
reference and the array output is: 

/N 

e^ = Z(q) - Z^(q) (A. 2) 

where e^ is a complex quantity. 

The mean square error is: 

®i®i = l®il^ " [Z(q) -Z^(q) J [Z(q) -Z^(q)J* (A. 3) 

le^l^ = lZ(q) 1^ -Z^(q)Z*(q) -Z*(q)Z(q) + |Z^(q)l^ (A. 4) 
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Substituting equation (A.l) into equation (A. 4) for Z^(q) 
yields ; 

|e^|^ = |Z(q) 1^ -Z*(q) (cTY(q)) -Z(q) (cTY(q))* + I A(q) 1^ 

(A. 5) 

The complex weight vector that will minimize equation (A. 5) 
is computed from the following recursive algorithm [Refs. 1, 
6,14] : 



= — i 2yj^ej^Y*(q) (A. 6) 

where is the convergence coefficient used in the gradient 
method [Ref . 6] . 

For the two-dimensional case discussed in Chapter II, 
Section 3.b, replace with cd^^ ^ to yield: 

“ cd + 2]jj^e^Y (q) (A. 7) 



where : 



e. = Z(q) - I I cd^ (m,n) Y (q,m,n) (A. 8) 

m n 



If the complex weights {cd^(m,n)} are separable, then 

cdj^(m,n) = c^(m)d^(n) 
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and : 



Ci_|_i (m) = c^(m) + 2y^e^(J (n) Y (q ,m, n) ) * (A. 9) 

n 

(m) = d^(n) + 2 (m) Y (q,m,n) ) * (A. 10) 

m 

where : 

e^ = 2(q) - I c^(m) I d^ (n) Y (q n) (A. 11) 

m n 

Reference 4 studied the application of the frequency domain 
LMS algorithm to a split array. The error signal is 
the difference between the outputs of the two arrays. It 
is possible to extend this idea to an (M xN) planar array 
by generating an error signal e^^(m,n) at each element location 
where : 



e^(m,n) = Z (q) - cd^ (m,n) Y (q,m,n) (A. 12) 

for all m,n. Although this scheme is not studied in this 
thesis, the implementation is rather straightforward. One 
needs only to replace the error e^^ with e^(m,n) in all the 
equations and to rewrite the code to generate the error . 

This scheme, however, does not have the advantage of 
averaging noisy data points from all of the elements in the 
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array to come up with an estimate. Therefore, it can be 
assumed that the noise performance will probably be inferior. 

The magnitude of the reference signal Z (q) is generated 
by averaging all (M xn) magnitude spectra whereas the phase 
of Z (q) is taken to be the phase of the reference element. 

A thresholding is done to determine which frequency bins 
should be processed. 

B. MODIFICATION TO THE LMS ADAPTIVE ALGORITHM TO MAKE IT 

USEFUL FOR SPATIAL RESOLUTION 

Three modifications are incorporated in the basic LMS 
adaptive algorithm. 

1 . Normalization of Recursive Complex Weights 

The components of the recursive complex weights (£j^/ 
and cd must be normalized to have magnitudes equal to 
one. This restriction forces the algorithm to achieve 
minimum error by adjusting the phase weights only. The esti- 
mation of direction cosines depends on the linear relationship 
of the phases of the signal across the array. Figure 35 
shows that, without the normalization performed after each 
complex weight update, the vector sum composing the estimate 
can be made close to that of the reference. The correct 
values of the direction cosines can only be achieved if each 
component in equation (A.l) has the same phase as the refer- 
ence phasor. All the vector sums shown in Figure 35 have 
the same resultant phase and magnitude as the reference but 
the phases of the components are not equal. It can also be 
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Figure 35. 



Different Ways of Adding Up Vectors 
to Obtain the Same Resultant 
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seen that there are infinitely many ways to achieve a small 
error without the normalization. The LMS adaptive filter 
implemented for this thesis has a unity magnitude across 
the array. It is, in fact, a phase-only filter. With the 
normalization, each component is forced to adapt to the 
phase angle of the reference. This is the only way that 
spatial resolution (based on the amount of phasor rotation 
of adaptive weights) can be achieved. 

2 . The Convergence Coefficient (y) 

This quantity is sometimes referred to as the feed- 
back coefficient, signifying the analogy to control systems 
[Ref,- 6]. For nonstationary environments, a constant con- 
vergence coefficient is usually chosen [Refs. 1,4,23]. 
However, if stationarity can be assumed, the performance of 
the LMS adaptive filter can be improved by using a mono- 
tonically decreasing scheme for implementing the convergence 
coefficient [Ref. 11]. If a constant convergence coeffi- 
cient is used, the bound on the choice of value is given 
by [Refs . 12,24]: 



0 < y < 



A 



-1 

max 



(A. 13) 



where 



A < 

max 



Tr [R] 



I I 

I R(i,i) = I R(0) 
i=0 i=0 



(I+1)R(0) 
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and R is the covariance matrix of the received noise cor- 
rupted signal. If y is larger than the maximum determined 
by equation (A. 13), oscillation will occur [Ref. 25] and 
the LMS adaptive filter will not converge. If y is too 
small, then the rate of convergence is slow. Figures 36, 37 
and 38 show the convergence characteristics for different 
values of y . 

The recursive LMS algorithm is based on the stochas- 
tic approximation technique [Ref. 11]. The choice of y is 
optimal if the following conditions are met [Ref. 12] : 



• y . > 0 

1 



lim y . = 0 

1 

l ->00 



I 



i=l 



r 2 

2 y. < “ 

i=l ^ 



A particular choice of y^ that meets the above four condi- 
tions is = 1/i [Ref. 12] . The effect of the adaptation 
decreases with the number of iterations and ceases completely 
for large i. Simulations show that the above choice of y^ 
is better than using a constant y, however, the choice 
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PERCENT ERROR IN ESTIMATING DIRECTION COSINE 
-120 -80 -40 



CONVERGENCE COEFFICIENT TOO LARGE 




Figure 36. Convergence Coefficient CConstant Value) 
Too Large 
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PERCENT ERROR IN ESTIMATING DIRECTION COSINE 
-80 -60 ’ -40 



CONVERGENCE COEFFICIENT TOO SMALL 




Figure 37. Convergence Coefficient (Constant 
Value) Too Small 
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PERCENT ERROR IN ESTIMATING DIRECTION COSINE 



BEST CONSTANT CONVERGENCE COEFFICIENT 




Figure 38. Optimum Convergence Coefficient 
(Constant Value) 
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(A. 14) 




where 0 < < 1 and is the initial value, appears to 

speed the rate of convergence and also achieve a high 

degree of accuracy. Equation (A. 14) does not meet the 

00 

condition J y. = ■», but simulations show that the response 
i=l ^ 

of the LMS adaptive filter is fast and accurate for 

0.75 < < 0.95. Further investigations of the convergence 

characteristics are warranted. The current implementation 

of the convergence coefficient results from experimentation 

with various values of a to achieve the best estimates of 

s 

direction cosines. With this scheme, the initial value of 
y^ can be set quite liberally since the value of y^ decreases 
geometrically . 

C. PHASE WRAP-AROUND PROBLEM 

The resolution of the phase wrap-around problem in the 
adaptive phase weights was discussed in Chapter II exten- 
sively. The scheme to resolve the phase ambiguity depends 
on the proper functioning of the elements adjacent to the 
reference element. In the event that some of those adjacent 
elements are not operational, the unwrapping scheme will 
have degraded performance. However, since the reference 
element can be any element in the array, it is possible to 
shift the location of the reference element to a region of 
the array that has sufficient adjacent elements that are 
functioning properly. 
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APPENDIX B 



DESCRIPTION OF SIMULATION PROGRAM 
FOR THE PASSIVE DETECTION CASE 



The VS APL application package worksapce (CMS file 

ADAPTIVE VSAPLWS) contains all the functions necessary to 

implement the simulation discussed in Chapter III. The 

general processing flow is as follows: 

- generate time samples of a plane wave signal of 

frequency f incident upon a (M xn) planar array with 
angles of incidence (0,(1)) 

add white Gaussian noise for a desired SNR 

compute the discrete Fourier transform (DFT) of each 
of the M xN time sequences 

determine the spectral line with the largest magnitude 
and its corresponding frequency bin number 

apply one of the three frequency domain LMS adaptive 
filters (orthogonal linear arrays, two-dimensional 
array with separable complex weights, or two-dimensional 
array) 

compute estimates of direction cosines. 

Usage of the functions are described below. 



A. SIG2D 

SIG2D generates planar array signal at each element 
location (equation (2.28)): 



Syntax: YN ^ AGl SIG2D AG2 



YN is a L xM xN matrix where L is the number of time 
samples, M is the number of elements in the x-direction, and 
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N is the number of elements in the y-direction. AGl is the 
spherical angle 0 and AG2 is the spherical angle 4> . The 
speed of sound (c ) , the signal amplitude (A) , the frequency 
of the signal (F ) , the interelement spacing (DX and DY ) , 
the sampling interval (T ) , the number of elements (M,N) 
and the number of time samples (L) can be changed by editing 
the function. 

B . NORRAND 

NORRAND generates independent white Gaussian noise 
samples . 



Syntax: NOISE ^ K NORRAND N Nl 

K is the number of noise samples desired, N is the mean 
Nl is the variance. The noise array NOISE must be reshaped 
to conform to the shape of the signal generated by SIG2G 
by : 



NOISE ^ L M N p NOISE 

The standard deviation = 1 is necessary to scale a sample func- 
tion of noise with zero mean and variance 1 to a desired 
signal-to-noise ratio. The signal power at each element is 
A /2 and the noise power at each element is The input 

signal-to-noise ratio at a single array element is then given 
by : 
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SNR 



(B.l) 





Solving for yields: 



2 1/2 



a 



N 



2 (SNR) ^ 



(B.2) 



where A is the amplitude of the signal and SNR is the numeri- 
cal signal-to-noise ratio. Therefore, a noisy signal with 
the desired SNR can be generated as: 



where NOISE has zero mean and variance 1. 

C . DFTWRT 

DFTWRT computes the discrete Fourier transform with 



YK has dimensions L xM xN; the first index L is nov/ the 
total number of frequency bins. RN is the noisy signal 
generated by adding noise to the output of SIG2D. A total 
of (M xN) L-point DFTs are computed using a radix-2 FFT 
algorithm. 



RN -f- YN + a., X noise 

N 



(B.3) 



respect to the time index for each element (equation (2.42)). 



Syntax ; YK DFTWRT RN 
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D . ADPLMS 

ADPLMS computes estimates of the direction cosines using 
the frequency domain LMS algorithm for the orthogonal array 
configuration . 



Syntax: QTGT ADPLMS YK 

QTGT is the frequency bin number where a valid signal has 
been identified and YK is the output of the function DFTWRT. 

A reference signal at frequency bin QTGT is generated by 
calling the function REFGEN. Direction cosine estimates 
u and V are computed in two different loops since in general 
the number of elements M is not necessarily equal to N. 
Estimates of both direction cosines are generated for each 
iteration under the names UHAT and VHAT. The phase unwrapping 
is done by calling the function DClDX for the x-direction 
and DClDY for the y-direction. The number of iterations 
(ITER) and the initial convergence coefficients (MUX, MUY) 
can be changed by editing the function ADPLMS. The scale 
factor (SCMU) for decreasing the convergence coefficients 
can be changed in the workspace. 

E . QFLMS 

QFLMS computes estimates of the direction cosines 
using the frequency domain LMS algorithm for the two-dimen- 
sional planar array with separable complex weights. 

Syntax: QTGT QFLMS YK 
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QTGT is the frequency bin number and YK is the output of the 
function DFTWRT. A reference signal is generated by the 
function REFGEN. Direction cosine estimates u and v are 
computed every iteration of the LMS adaptive loop. They 
are stored in the vectors named UHATQF and VHATQF. The 
phase unwrapping is done by using DCIDX and DClDY in the x 
and y-directions , respectively. Recall that only M+N com- 
plex weights are updated for this configuration because of 
the separability assumption. The initial convergence 
coefficient (MUQ) and the number of iterations (ITER) can 
be changed by editing the function QFLMS. The scale factor 
is named SCMU and is stored in the variable list in the 
workspace. 

F. ADPLMS2D 

ADPLMS2D computes estimates of the direction cosines 
using the frequency domain LMS adaptive algorithm for a 
two-dimensional planar array. 

Syntax: QTGT ADPLMS2D YK 

QTGT and YK are the same quantities described in the last 
two functions. A reference signal at frequency bin QTGT 
is generated. Direction cosine estimates u and v are com- 
puted for each iteration and stored in vectors named UHAT2D 
and VHAT2D. The phase unwrapping is accomplished by using 
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the function DC2D. The number of iterations (ITER) and the 
initial convergence coefficient (MUG) can be changed by 
editing the function ADPLMS2D. The scale factor can be 
changed by assigning a different value to SCMU in the work- 
space. The use of a different unwrapping function is 
required since the complex weights are not assumed to be 
separable. 

G . REFGEN 

REFGEN generates a reference signal at a particular 
frequency bin QTGT. 

*■' • 

Syntax: CQREF ^ QTGT REFGEN YK 

CQREF is the reference signal used in the frequency domain 
LMS adaptive filter. The magnitude of CQREF is obtained by 
averaging the magnitudes of all (M xn) frequency spectra 
in the frequency bin QTGT. The phase of CQREF is taken 
to be the phase of the reference element in the QTGT fre- 
quency bin. 

H. DClDX 

DClDX unwraps the phase weights for a linear array in 
the x-direction. 

Syntax: UHAT[i] ^ N DClDY DV 
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UHAT[i] is the direction cosine estimate u of the i 
iteration, M is the number of elements in the x-direction 
and CV is a complex vector of M phase weights that will co- 
phase the incident signal. 

I. DCIDY 

DCIDY unwraps the phase weights for a linear array in 
the y-direction. 

Syntax: VHAT[i] N DCIDY DV 

th 

VHAT[iJ is the direction cosine estimate v of the i iter- 
ation, N is the number of elements in the y-direction and 
CV is a complex vector of N phase weights that will cophase 
the incident signal.. 

J. DC2D 

DC2D unwraps the phase weights for a two-dimensional 
array . 



Syntax: DC2D CD 

CD is the two-dimensional phase matrix that will cophase 
the incident plane wave signal. This function is used by 
the function ADPLMS2D. 

An example of using this application package is shown 
as follows: 
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YN ^ 



55 SIG2D 35 



signal generation 



NOISE ^ 3200 NORRAND 0 1 



noise generation 



NOISE 128 5 5 p NOISE 



noise generation 



RN -s- YN + SCALE NOISE 



noisy signal 



YK ^ DFTWRT RN 



transform to frequency 
domain 



QTGT ADPLMS YK 



estimate direction 
cosines 



QTGT ADPLMS 2d YK 



estimate direction 
cosines 



QTGT QFLMS YK 



estimate direction 
cosines 



This sequence of statements generates a plane wave signal 
incident upon a 5 x 5 planar array with angles of incidence 
6 = 55° and (p = 35°. The number of time samples for each 
element is 128. A noise matrix is then generated and added 
to the plane wave signal and the discrete Fourier transform 
with respect to time is taken for each element in the array. 
The three array configurations for the frequency domain LMS 
algorithm are then used sequentially to estimate the 
direction cosines of the incident plane wave. 
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APPENDIX C 



DESCRIPTION OF SIMULATION PROGRAM 
FOR THE PULSE COMMUNICATION CASE 

The VS FORTRAN prgrams were written to implement the 
pulse communication problem discussed in Chapter IV. Two 
separate programs were written utilizing essentially the 
same subprograms. These programs handle the quadrature 
demodulated complex envelope signals generated by Vos' [Ref. 
19] program. The programs are available on user account 
0218P at the Naval Postgraduate School, Monterey, California. 

A. PROGRAM ADBFM 

This program is compiled using FORTVS and is designed 
to run under DISSPLA. It requires a storage capacity of 
1 Mbyte. The following sequence of commands should be used 
to run the program. 

- DEFINE STORAGE 1 M 

- I CMS 

- GLOBAL TXTLIB VALTLIB VFORTLIB CMSLIB IMSLSP NONIMSL 

- GLOBAL LOADLIB VFLODLIB 

- FILEDEF 04 DISK fname DATA 

(fname is the filename of the date file) 

- LOAD ADBFM 

- DISSPLA ADBFM 

When execution begins, the user will be prompted to enter 
the desired values for the necessary parameters. These 
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parameters are noise status, input signal-to-noise ratio in 
dB at a single array element, number of iterations, spec- 
tral line to be processed, convergence coefficient, scale 
factor for the convergence coefficient, and the choice of 
one of the three array configurations. 

For each array configuration, plots of the estimates of 
the direction cosines are generated. The plots for magni- 
tude and phase of the difference between the reference 
signal and the estimate are also generated. 

B. PROGRAM ERVSDB 

This program computes the rms errors for various input 
signal-to-noise ratios. If a plot of rms error versus SNR 
(dB) is not required, this program does not have to be run 
under DISSPLA. The following sequence should be used: 

- DEFINE STORAGE 1 M (if plot is required) 

- I CMS 

- GLOBAL TXTLIB VALTLIB VFORTLIB CMSLIB IMSLSP NONIMSL 

- FILEDEF 04 DISK fname DATA 

(fname is the filename of the data file) 

- LOAD ERVSDB 

- DISSPLA ERVSDB (if plot is needed) 
or START * (if plot is not needed) 

When execution begins, the program will prompt the user to 

enter an initial input signal-to-noise ratio in dB at a single 

array element. It will then ask for a dB step size such 

that the next SNR is determined by the current SNR in dB 

plus the dB step size. A total of nine SNRs are allowed. 
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For example, if the initial SNR is -12 dB and the step size 
is 3 dB, then the program will compute rms errors for the 
set of dB levels (-12, -9, -6, -3, 0, 3, 6, 9, 12}. The 
other parameters such as iteration number, initial conver- 
gence coefficient, scale factor for convergence coefficient, 
and the spectral line to be processed are entered when 
prompted. The program will then ask for how many sample 
functions of signal and noise are to be averaged. Simula- 
tion results show that the average of 50 to 100 sample func- 
tions are sufficient to reduce the variance of the direction 
cosine estimates. One of the three array configurations is 
then -chosen by the user to estimate the direction cosines. 

The screen output of this program is ordered pairs of 
rms errors corresponding to a particular input SNR in dB 
for all nine specified SNRs . A plot of rms error versus 
SNR in dB can be generated if desired (provided that the 
program is run under DISSPLA) . 
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